function [out,outcell,hand_vec,sumSt]=dsge_spectrum(GG,RR,ZZ,SDX,specrep,stanames,shonames,serdesc,lowper,highper,npoints, outpath, figname) 
%
% [out,outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,specrep,stanames,shonames,serdesc,lowper,highper,npoints,outpath) 
%
% =========================================================================
% DSGE_SPECTRUM 
% Compute the Spectral Decomposition of a DSGE, reporting selected states 
% raw or in levels, as well as the observables. 
% 
% SPECREP   Names of states wish to report. If ' lev' then it will be 
%           multiplied by the (1-L)^(-1) filter 
% STANAMES  Cell with the names of all states 
% SHONAMES  Cell with the names of all shocks 
% SERDESC   Cell with the names of observables, in the ZZ matrix 
% LOWPER    Lower periodicity > 2 
% HIGHPER   Higher periodicity 
%           These two determine the frequency band over which the spectrum
%           is computed. For business cycles, usually [6 32] are used        
% NPOINTS   Number of points used in the approximation, 1000 should provide
%           a reasonable approximation 
% OUTPATH   If provided, a table and figures will be stored 
% 
% Alejandro Justiniano April 28 2010 
% ========================================================================

% 1. OBTAIN Indicators of which rows of the state vector to report and
% whether to transform to a level or not 
% =========================================================================
% STSPEC        Position in STANAMES of the states to report 
% SPECLEVIND    [0 1] Indicator of which of these states will be in levels
% STNAMERAW     Is the name of the RAW, i.e. untransformed states 
if nargin < 12 || isempty(outpath); 
    flag_save=0; 
    hand_vec=[]; 
else 
    flag_save=1; 
end 

stanames=stanames(:); 
shonames=shonames(:); 

nspec     =length(specrep);
stspec    =zeros(nspec,1);
speclevind=zeros(nspec,1);
stnameraw =cell(nspec,1); 

NLevel=0; 
for ii=1:nspec;
    tempos=findstr(' lev',specrep{ii} );
    if ~isempty( tempos);
        tempstr=specrep{ii};
        stnameraw(ii)={tempstr(1:tempos-1)}; 
        stspec(ii) =cellposition( stnameraw(ii) , stanames );
        speclevind(ii)=1;
        NLevel=NLevel+1; 
    else
        stspec(ii) =cellposition( specrep(ii), stanames );
        speclevind(ii)=0;
        stnameraw(ii)=specrep(ii); 
    end
end
clear tempos tempstr;
poslevind=find( speclevind == 1); 
disp('Variables for which to compute a Spectrum'); 
printcell( [stanames(stspec) num2cprec(speclevind)] ) 

% 2. Create frequency vector and associated bin-withds 
%    Note: the size of the frequency vector is NPOINTS 
%    If wish to use only the harmonic frequencies, according to the sample
%    size, use FREQ_VEC instead and let NPOINTS=size(Y,1) 
% =========================================================================
% a. Vector of peridicities 
pervec  =linspace(lowper,highper,npoints+2);
% b. Vector of frequencies 
freq    =fliplr((2*pi)./pervec);
binvec  =freq(2:end)-freq(1:end-1);
freq    =freq(2:end-1);

pervec  =(2*pi)./freq; 
% Obtain distance to edges 
% FREQ is in ascending order 
% WLB and WUB are the upper and lower bound of the frequencies 
% which are guaranteed to be within the desired set 
% EDGEW adjusts the rectangles at the end by the distance to the edges 
wlb  =2*pi/highper;wub=2*pi/lowper; 
edgew=[(freq(1)-wlb) (wub-freq(end))];
if any( edgew < 0 ); error('EDGEW cannot be negative'); end 

% Harmonic Frequencies 
% [freq,perv,edgew]=freq_vec(lowper,highper,npoints)

% COMPUTE the OUT structure 
% -------------------------------------------------------------------------
disp('Begin Spectral Analysis....');
dispaj('For periodicities ',['(',num2str(lowper),',',num2str(highper),')']); 
out=vardecom_spec(GG,RR,SDX,ZZ,freq,edgew,stspec,speclevind); 
%% Make Plots for the Variance Decomposition at different frequencies - Leonardo 10/8/2018
%makeplots
disp('Finished Spectral Analysis' ); 

% 4. SPECDEC    Variance decompsition for that frequency
%               Diagonal elements ONLY
% -------------------------------------------------------
% SPECDECZ      [NZ NX NF]
% SPECDECST     [NST NX NF]
% SPECDECLEV    [NLEV NX NF]
%

% OUTPUT 
% OUT reports 
% 4. SPECDEC    Variance decompsition for that frequency
%               Diagonal elements ONLY
% -------------------------------------------------------
% SPECDECZ      [NZ NX NF]
% SPECDECST     [NST NX NF]
% SPECDECLEV    [NLEV NX NF]
%

% I. Plots showing the variance shares attributed to each shock 
%    plotted versus the periodicities 
% 
% II. Tables showing the average decomposition over the frequency band 
% 
% Report 
% 5.VFBDECL total Variance decompositon for this frequency band
% Ratio of the total variance in SPECDEC to the total variance in
% VFB
% ------------------------------------------------------------
% VFBDECZ      [NZ NX]
% VFBDECST     [NST NX]
% VFBDECLEV    [NLEV NX]


% I. Top cell with settings followed by decompositions 
% =========================================================================
empstr=['(',num2str(lowper),',',num2str(highper),')']; 
outcell=emptycell(2,6); 
outcell(1,[1:2  4:5])={'Lower Per',num2str(lowper),'npoints',num2str(npoints)}; 
outcell(2, 1:2     )={'High  Per',num2str(highper)}; 
NSummary=length(stnameraw)+length(serdesc)+NLevel; 


sumSt.mat=zeros(NSummary,length(shonames)); 
sumSt.names=cell(NSummary,1); 

if ~isempty(poslevind)
    sumSt.mat(1:NLevel,:)=out.vfbdeclev; 
    sumSt.names(1:NLevel)=specrep(poslevind); 
    summaryPos=NLevel+1; 
    outspec=crtcell( out.vfbdeclev ,  specrep(poslevind) , shonames ,{'Spec Dec LEVELS'} );
    outcell=merge_cells( outcell, outspec , 1 );
else 
    summaryPos=1;
end

outspec=crtcell( out.vfbdecst  ,  stnameraw          , shonames ,{'Spec Dec ST Select'} );
rowPos=summaryPos:(summaryPos+length(stnameraw)-1);
sumSt.mat(rowPos,:)=out.vfbdecst;
sumSt.names(rowPos)=stnameraw;
summaryPos=rowPos(end)+1;

outcell=merge_cells( outcell, outspec , 1 ); 
outspec=crtcell( out.vfbdecz   ,  serdesc            , shonames ,{'Spec Dec Z        '} ); 
rowPos=summaryPos:(summaryPos+length(serdesc)-1);
sumSt.mat(rowPos,:)=out.vfbdecz;
sumSt.names(rowPos)=serdesc;

outcell=merge_cells( outcell, outspec , 1 ); 
if flag_save==0 
    return 
end 
% II. Plots of decompositions for the states 
%%
nrow=ceil( length(shonames)/2 ); 
ncol=2; 
countlev=0; 
pervecflipped=fliplr(pervec); 
xaxis=round(linspace(floor(pervecflipped(1)),ceil(pervecflipped(end)),5));

hand_vec=zeros(nspec,1); 

titfont=12; 
for ii=1:nspec; 
    hand_vec(ii)=figure;    
    if speclevind(ii)==1
        countlev=countlev+1;
    end
    
    for jj=1:length(shonames);
        subplot(nrow,ncol,jj);
        
        if speclevind(ii)==1;
            plot( pervecflipped, flipud( squeeze( out.specdeclev(countlev,jj,:) ) ), '-b','LineWidth',2.5 );
            yvec=zeros(2,1); 
            yvec(1)=max( min( squeeze( out.specdeclev(countlev,jj,:) )  -0.01 ), 0 ); 
            yvec(2)=min( max( squeeze( out.specdeclev(countlev,jj,:) )  +0.01 ), 1 ); 
        else
            plot( pervecflipped, flipud( squeeze( out.specdecst(ii,jj,:) ) ), '-b','LineWidth',2.5 );
            yvec(1)=max( min( squeeze( out.specdecst(ii,jj,:) )  -0.01 ), 0 );
            yvec(2)=min( max( squeeze( out.specdecst(ii,jj,:) )  +0.01 ), 1 );
        end
        ylim([yvec(1) yvec(2)]); 
        xlim([pervecflipped(1)-1 pervecflipped(end)+1]);         
        set(gca,'Xtick',xaxis);
        title(shonames{jj},'FontSize',titfont);
        
        if jj==length(shonames); 
            xlabel('Periodicity'); 
    ylabel('Share      '); 
   
        end 
    end
    
    suptitle( [upper(specrep{ii}),':  Spectral Decomposition  ',empstr], 14 );     
    set(hand_vec(ii),'PaperPosition',[0.75 0.5 7.5 10])
    
    subtitulo( outpath); 
    
end 
%%

cucd=pwd; 
if nargin < 13; 
    figname = 'spectral decomposition fig'; 
end 

if flag_save==1
    cd(outpath) 
    print(hand_vec(1),'-dpsc',figname) 
    if nspec > 1
        for ii=2:nspec         
            print(hand_vec(ii),'-dpsc',figname,'-append') 
        end 
    end
    if isunix == 0
        xlswrite('variance decomposition tab',outcell,'spectrum'); 
    else
        savecell(outcell,[],outpath,   'variance decomposition tab');   
    end
     
    if isunix == 0
        ps2pdf('psfile', [figname,'.ps'], 'pdffile', [figname, '.pdf'])
    end
end 
cd(cucd) 
